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The development and verification of the Charring Ablating Thermal Protection Implicit System Solver is 
presented. This work concentrates on the derivation and verification of the stationary grid terms in the equa- 
tions that govern three-dimensional heat and mass transfer for charring thermal protection systems including 
pyrolysis gas flow through the porous char layer. The governing equations are discretized according to the 
Galerkin finite element method with first and second order implicit time integrators. The governing equations 
are fully coupled and are solved in parallel via Newton’s method, while the fully implicit linear system is solved 
with the Generalized Minimal Residual method. Verification results from exact solutions and the Method of 
Manufactured Solutions are presented to show spatial and temporal orders of accuracy as well as nonlinear 
convergence rates. 


Nomenclature 


a thermal diffusivity ( m2 /sec) or absorptivity 

a t , f3 t , 7 t temporal finite difference weights (sec -1 ) 

/ 3 extent of reaction 

Q volumetric flow rate ( m3 /sec) 

U independent variable vector 

v velocity ( m / sec) 

v' superficial velocity ( m / sec) 

Ah characteristic mesh length scale (m) 

At time step (sec) 

s surface recession rate ( m /sec) 

6 j mass source ( k s/ m 3 -sec) 

171 mass flux ( k g/m 1 2 -sec) 

<5 energy source ( w /m 3 ) 

q heat flux ( w /m 2 ) 

e emissivity or penalty factor 

r volume fraction in virgin material 

h unit normal vector 

k, k. permeability (m 2 ) 

1Z residual 

p dynamic viscosity (Pa • sec) 

</> porosity 

if) basis function 

p density ( k s/ m 3 ) 

er Stefan-Boltzmann constant ( w /m 2 K 4 ) 

6 non-dimensional temperature or manufactured solution amplitude 
A area (m 2 ) 

A tl B t ,Ct, D t manufactured solution temporal coefficients (sec -1 ) 
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A x , By, C z manufactured solution spatial coefficients (m -1 ) 

Ch heat transfer Stanton number 

C p specific heat ( J /kg K) 

DE discretization error 

E activation energy ( J /kg) or error 

e internal energy ( J /kg) 

e G total internal energy ( J /kg) 

H heat transfer coefficient ( w /m 2 K) 

h enthalpy ( J /kg) 

h° enthalpy of formation ( J /kg) 

h a total enthalpy ( J /kg) 

J x W Jacobian of the transformation matrix times the integration weight 
k. k pre-exponential factor (sec -1 ) or thermal conductivity ( w /m K) 

L domain length scale (m) 

to reaction order 

N number of nodes 

N e number of elements 

nc number of components 

P pressure (Pa) 

p temporal order of accuracy 

Q generic quantity of interest 

q spatial order of accuracy 

R gas constant (J/kg-K) 

r perturbation factor 

T temperature (K) 

t time (sec) 

u velocity component ( m /sec) 

v test function 

Subscripts and Superscripts 
oo denotes far-held state 

v iteration index 

c quantity of the char 

e boundary layer edge or element quantity 
g quantity of the gas 

h denotes finite dimensional approximation 

to quantity of the mesh 

n time index 

o denotes total, initial, or reference condition 
QP denotes a value at a quadrature point 
r denotes recovery state 

s quantity of the solid or surface 

v quantity of the virgin 

w quantity of gases adjacent to the wall 


I. Introduction 

Reliable predictions of the performance of thermal protection systems (TPS) is critical in the design and analysis 
of atmospheric entry vehicles. For decades, design engineers have relied on numerical heat transfer solutions for 
the highest fidelity thermal response predictions. For vehicle components exposed to high heat fluxes, greater than 
40 w /cm 2 , usually an ablative thermal protection system is required, and they generally tend to be more thermally and 
chemically complex than reusable TPS materials. These complexities require additional physics to be accounted for 
in the computational models. 

In general, there are two classes of ablators: those that only undergo surface ablation and those that decompose, 
char, or pyrolyze internally in addition to ablating at the surface, which is depicted in Figure 1 . Before entering an 
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atmosphere, a charring ablator TPS material initially starts out as the virgin composite material. The virgin composite 
is typically a fibrous material injected with resin that binds the material and adds mechanical strength. As the vehicle 
enters, it is exposed to a heat flux which leads to thermal penetration in-depth. As the material heats up, resin and 
fiber constituents can pyrolyze, which is typically an endothermic process, resulting in voids and gas generation. The 
pyrolysis gases are driven through the porous char layer by internal pressure gradients, and the majority of the gas 
is expelled from the surface into the chemically reacting boundary layer. As the surface temperature rises, boundary 
layer gases start to react with the exposed char and pyrolysis gas at the surface. The gas/surface chemical reactions 
result in mass removal from the surface and the ablation products are carried away with the fluid stream. This mass 
removal causes the surface to recede, and the ablative TPS thins throughout entry. 



Figure 1. Charring ablator schematic. 


There have been many attempts to model the complicated thermochemical processes involved in the charring 
ablator problem for the past fifty years. The most notable contributions to this technology area were by the Aerotherm 
Corporation in the late 1960s. They presented the development of the charring material, surface thermochemistry, and 
reacting boundary layer models used in CMA, ACE, and BLIMP respectively. 1 5 The CMA theory manual written by 
Moyer and Rindal is of particular interest to this study since it outlines the model for the one-dimensional energy and 
mass transfer problem within the ablator. Moyer and Rindal proposed an in-depth material model that is still being 
used today. The model implemented in CMA assumes that the virgin and char thermal properties are known, and the 
thermal properties in the decomposing region is interpolated between the two. The decomposition of virgin into char 
and gas is modeled using Arrhenius rate equations. The pyrolysis gases resulting from decompostion of the virgin 
material have no residence time in the char layer and are instantaneously injected into the boundary layer as soon as 
they are produced. The pyrolysis gases do not react with the char layer as they percolate to the surface nor do they 
deposit residue (coke) in the pore space. 

Numerous researchers and code developers have adopted the material model in CMA since its introduction and 
have made incremental improvements including improved mesh motion schemes over the translating/node-dropping 
scheme in CMA, improved discretization methods and numerical algorithms, extension to multiple dimensions, ad- 
dition of gas phase continuity equation and porous flow equations, and improved physical models such as thermo- 
chemical non-equilibrium in-depth, del Valle, Pike, and April developed a model in which they solved the equilibrium 
and non-equilibrium porous gas flow equations to determine the pressure and gas composition in-depth. 4 5 Chen and 
Milos developed the FIAT code, which is one dimensional and has the same material model as CMA, but it is fully im- 
plicit. 6 FIAT has been used in numerous TPS design and analysis activities since its development which most recently 
includes Stardust, Mars Science Laboratory, and Orion. 7 9 Blackwell, Hogan, and Cochran introduced the use of the 
control volume finite element method (CVFEM) and linear elastic mesh motion to solve one and two dimensional 
ablation problems. 10-12 Ahn and Suzuki et. al. developed the SCMA and SCMA2 codes that solve the one and two 
dimensional charring ablation problems including pyrolysis gas flow within the char layer. 1314 Amar, Blackwell, and 
Edwards used the CVFEM to solve the coupled energy and gas continuity equations and presented verification studies 
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which proved correct implementation of many of the boundary conditions and terms in the governing equations. 15-17 
Recently, Dec and Braun presented codes that solve one and three-dimensional ablation problems including pyrolysis 
gas flow. Dec and Braun used the Galerkin finite element method (FEM) for spatial discretization and use linear elastic 
mesh motion. 18,19 

One recurring shortcoming of most of the ablation code development work over the years has been the lack of code 
verification studies to ensure proper implementation of the governing equations, boundary conditions, and numerical 
methods employed in the software. While various definitions for “code verification and validation” can be found in 
literature, the definition for code verification adopted in this work is given by Roache: 20 

The code author defines precisely what continuum partial differential equations and continuum bound- 
ary conditions are being solved, and convincingly demonstrates that they are solved correctly, i.e. usually 
with some order of accuracy, and always consistently, so that as some measure of discretization fe.g. 
mesh increments) A — > 0, the code produces a solution to the continuum equations; this is Verification. 
Whether or not those equations and that solution bear any relation to a physical problem of interest to the 
code user is the subject of Validation. 

In other words, verification is the strictly mathematical exercise of proving that the model equations and numerical 
methods have been correctly implemented in the software, and validation provides evidence that the correct equa- 
tions have been chosen to solve the physical problem of interest. While many code developers bypass verification 
exercises and attempt to validate their codes through comparisons with data, verification is a necessary step that must 
precede code validation because one can not know if the correct equations are being solved (validation) without first 
knowing that the equations they intend to solve are correctly implemented (verification). Additionally, code-to-code 
comparisons and/or extensive use of a code are not part of the rigorous verification process. However, these activities 
do have merit in that they help build confidence, improve code usability, and demonstrate robustness and reliability. 
While verification studies should be performed for every term in the governing equations, all boundary conditions, 
and all modeling options, there is some practical limit to how rigorously a developer can verify a code given that most 
software of this nature is developed in an environment where sooner than later practical problems need to be solved 
for engineering design and analysis. 

In this study, the development and verification of the Charring Ablating Thermal Protection Implicit System Solver 
is presented. The code solves the two and three dimensional charring ablator, general porous flow, and heat transfer 
problems in serial or parallel. The code uses the first and second order implicit time integrators and the Galerkin FEM 
spatial discretization on unstructured hybrid meshes. Supported element types include quadrilaterals and triangles in 
two dimensions and tetrahedrons, hexahedrons, prisms, and pyramids in three dimensions. Furthermore, the code is 
capable of solving problems using adaptive mesh refinement on all element types in both two and three dimensions. 
The governing equations are fully coupled and solved fully implicitly with Newton’s method (with exact Jacobian 
terms), and the Generalized Minimal Residual method (GMRES) is used to solve the linear system. 21 The exact 
Jacobian terms are calculated through the use of the complex-step (or complex perturbation) method rather than being 
analytically derived. This enables optimal convergence of the nonlinear solver as opposed to degraded convergence 
rates from approximate Jacobians. As far as the authors are aware, the complex-step method is a novel approach to 
calculating Jacobians for numerical heat transfer or CFD. The method dramatically reduces development effort and 
the potential for bugs resulting from derivation and coding errors. 

This study concentrates on the verification of the energy and porous flow equations in both two and three dimen- 
sions, several of their boundary conditions, and the nonlinear system solver. Comparisons with analytical solutions 
are used when available but the Method of Manufactured Solutions (MMS) is also used to verify the nonlinear heat 
and mass transfer equations. 20,22 This study only deals with the non-decomposing and stationary mesh terms. While 
not all of the physical models used in this study are necessarily relevant for a broad range of ablation and porous flow 
problems, it should be realized that the terms in the governing equations verified here will be present in all charring 
ablation problems. Consequently, these are the necessary first steps in having a fully, or at least adequately, verified 
code. 


II. Mathematical Model 


II.A. Governing Equations 

The equations that govern the solid/gas system of the porous charring ablator include energy and mass conservation 
equations for the solid as well as the Navier-Stokes equations as applied to all of the gaseous species considered. In 
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the general case, it is possible that the pyrolysis gases react with the remaining solid, or deposit residue (coke) on 
the solid, but these phenomena are neglected. Under the assumptions that the pyrolysis gas is in thermochemical 
equilibrium and the solid and gas are in thermal equilibrium, then the solid and gas energy equations for a moving 
mesh reduce to a mixture energy equation given by 

= v • (fevr) - V • (</> PghoVg ) + v • (ph 0 v m ) + Q (1) 


The velocity of the pyrolysis gases is governed by a porous flow law such as Darcy’s law, which omits the necessity 
to solve the momentum equations. Since ablators in general can be anisotropic materials, the thermal conductivity is 
a second order tensor. The solid mass conservation equation is 

^ = Cj s + V- ( p s V m ) (2) 

If it is assumed that all solid decomposition results in pyrolysis gas generation, aj g = —u) s , and the gases occupy all 
of the pore space and are free to flow through it, then the gas mass conservation equation is given by 

= w a - V • {<j>p g v g ) + V • {<f>p a v m ) (3) 

In this study, not every term in the governing equations will be verified. Looking only at the terms that are verified in 
this work, the equations essentially describe the energy and mass transfer in a non-ablating material (no mesh motion) 
assuming no in-depth decomposition, no gas flow contributions in the energy equation, and no energy contribution in 
the gas flow equation. Consequently, the verified equations reduce to 


and 


Verified Energy Equation: 


= v • (fevr) + q 


(4) 


Verified Gas Continuity Equation: 


9(<PPg) 

dt 


V • {(jtpgVg) 


(5) 


Although only a subset of the terms in the governing equations are verified in the current work, every term has been 
implemented in the code. Consequently, the material model and boundary conditions used in the code will be presented 
even though not all parts of the model are exercised in the problems presented here. 


II.B. Physical Boundary Conditions 

For most entry type ablation analyses, the incident heat flux is typically a combination of convective and radiative 
components. If a film coefficient boundary layer model is adopted where the species diffusion coefficients are assumed 
equal, and the mass and heat transfer film coefficients are assumed equal, then the energy balance for an infinitesimally 
thin control volume attached to the receding surface reduces to 

Convection Radiation Reradiation Wall Flux Char Flux Gas Flux Conduction 

Pe^eC H {hr h u i) T txq ra d (76 {T. w ^oo) ht w h 0w lfl c h c T tflg s h 0g T Qcond s — 0 (6) 

and the surface mass balance is 

rhg s — m w — fn c = 0 (7) 

where 

Wall mass flux injected into fluid: rh w = p w v w ■ ft (8) 

Solid (char) flux to the surface from surface recession: fn c = p c s ■ ft (9) 

Pyrolysis gas flux to surface from in-depth convection and surface recession: m gs = (j>p g {v g — s) ■ ft (10) 

Conduction flux on solid side: q con d s - —k s 'VT s ft (11) 

and n is a unit normal vector pointing outward into the fluid. The boundary condition for the gas continuity equation 
is typically a specified pressure condition. Since there are no moving mesh ablation problems in the current study, 
several other simpler boundary conditions are verified. For the energy equation these include: 

Specified input heat flux: q spe c + qcond s =0 (12) 
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Perfect gas convection: H ( T r — T w ) + q con d s =0 (13) 

and specified temperature. For the gas continuity equation, the boundary conditions verified in this study are specified 
pressure, specified density, and 


Specified inflow mass flux: ih gspec + m gs = 0 (14) 

Since there is no surface recession, the surface gas flux in Eq. (10) becomes 

rn.g s = (\>pgV g -h (15) 

The implementation of these boundary conditions in the finite element formulation will presented in a later section. 


II.C. Material Model 

In order to sufficiently explain the governing equations and boundary conditions, it is important to understand the ma- 
terial model used to characterize the state of the solid/gas mixture. It is assumed that all the pores are interconnected, 
and therefore the gas occupies all of the pore space and is free to flow through it. Consequently, the density of the 
solid/gas mixture is described by 

p = (j)fj g + p s (16) 

where the solid density is a bulk density, the gas density is a density with respect to the space the gas occupies (pore 
space), and the porosity is equal to the gas volume fraction. In terms of units, Eq. (16) can be expressed as 


p 

^ A N 

[total mass] 
[total vol] 


<t> Pg 

[pore vol] [gas mass] 
[total vol] [pore vol] 


Ps 

[solid mass] 
[total vol] 


(17) 


It is assumed that the thermodynamic state of the pyrolysis gases can be described as a mixture of perfect gases, and 
that the solid and gas phases are in thermal equilibrium resulting in 


T g = T s =T 


(18) 


P = f(p g ,T ) (19) 

where the equation of state and all relationships describing the thermochemical state of the gas are dictated by an 
external equilibrium chemistry library such as Cantera 23 or CEA 24 ' 25 to which the code is linked. During the code 
development effort, Cantera has been used to provide the thermochemical state of the gas. It should be noted that 
assuming that the gas is in thermochemical equilibrium may be a poor assumption for some problems. However the 
equilibrium assumption is being used for development purposes. Using a library such as Cantera does allow for easy 
extension to chemical non-equilibrium modeling. 

The solid material model adopted in this study is similar to the model developed by Moyer and Rindal 1 but has 
been expanded to include an arbitrary number of components, nc. The solid bulk density is given by 


= TlPi ( 20 ) 
i=l 


where T,; is the volume fraction of the i th component in the virgin composite and is therefore constant. The units 
associated with the solid bulk density model in Eq. (20) are 


Ps 

, * V 

[solid mass] 


[total vol] 


E 


r* P i 

, ^ ^ ^ , 

[initial vol of i th comp.] [mass of i th comp.] 

[total vol] [initial vol of i th comp.] 


( 21 ) 


It is assumed that the solid does not change volume due to thermal expansion, and therefore the total volume is 
constant. It is important to note that the solid description in Eq. (20) is only a modeling assumption, and the solid is 
not truly comprised of nc components, species, or distinguishable materials. This modeling assumption comes as a 
result of decomposition data obtained from thermogravimetric analyses (TGA). 
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It is assumed that all decomposed solid mass results in gas generation, and the general model of the decomposition 
process is described by the irreversible reaction 


virgin 


char + gas 


Taking the temporal derivative of Eq. (20) gives the solid decomposition rate in terms of component decomposition 
rates. 


A TIC A 

dPs _ ^ p OPi 


at 


dt 


( 22 ) 


It is assumed that the decomposition of each component can be described by an Arrhenius relationship of the form 


dpi 

dt 


kiPvi 



rrii 

e -Ei/RT 


(23) 


which applies at a constant spatial location (as apposed to a given node which can move during the solution process). 

Since most thermophysical properties of the solid are only known for the virgin and fully charred states, the inter- 
mediate solid is modeled as some interpolated state between virgin and char. This interpolated state is characterized 
by the extent of reaction (T), or degree of char, given by 


0 = (24) 

Pv pc 

where the virgin and char bulk densities are known constants. It is evident that as the solid decomposes from virgin to 
char, the extent of reaction ranges from 0 to 1. The definition in Eq. (24) can be rearranged to more clearly describe 
the interpolated state. 

Pa = (1 - P) Pv + Pp c (25) 

Although the virgin and char materials are not distinguishable entities within the intermediate solid, Eq. (25) reveals 
that the degree of char represents an effective char volume fraction within the solid (not in the solid/gas mixture). 


II.D. Property Models 

II. D. 1. Internal Energy and Enthalpy 

The total internal energy and enthalpy of the system can be described by 


pe Q — (1 f3) Pv^-v T dp,fic T tf)pg€ 0g 


(26) 


and 

where 


ph a (1 4 ) PiPv 4 ~ PPcPc 4 ~ (/>PgK g 

e v /c = h v/c = h° v/c + [ C Pv/c (T')dT' 
Jt 0 


and the total energy of the gas is 


Similarly, the total enthalpy of the gas is 


1 / 

e O g =e g +-{Vg- Vg) 

K g =hg+^(Vg- Vg) 


(27) 

(28) 

(29) 

(30) 


where the specific internal energy and specific enthalpy come from the chemistry library. Notice that the kinetic energy 
of the gas is not neglected in the current model in contrast to the CMA model. While it may be appropriate to neglect 
the kinetic energy for many ablation problems, it may not be appropriate for other porous flow problems that this code 
is intended to solve. 
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II.D.2. Thermal Conductivity 


The thermal conductivity is in general an anisotropic function of temperature and pressure resulting in a dim-by- 
dim tensor where each component of the tensor can be an independent function of temperature. For materials with 
anisotropies not aligned with coordinate axes, an additional transformation matrix must be introduced to appropriately 
apply the model. The thermal conductivity of the mixture is assumed to be a mass weighted average of virgin char and 
gas. 


k = 


(1 - /?) p v - 

K v 

p 




(31) 


While more appropriate combination laws may exist to determine the mixture thermal conductivity, this simplified 
model has been chosen for development purposes. 


II. D. 3. Permeability 

The permeability is in general an anisotropic function of the extent of reaction defined by Eq. (24) resulting in a 
dim-by-dim tensor where each component of the tensor can be an independent function of extent of reaction. 

k = k{(3) (32) 

Like the thermal conductivity, for materials with anisotropies not aligned with coordinate axes, an additional transfor- 
mation matrix must be introduced to appropriately apply the model. 

II.D.4. Porosity 

The porosity is in general a function of extent of reaction, and represents the gas mass fraction since all of the pore 
space is assumed to be interconnected. 

<!> = HP) (33) 


II. D. 5. Emissivity and Absorptivity 

The emissivity and absorptivity of the mixture are assumed to be equal to the emissivity and absorptivity of the solid. 
They are modeled similarly to the thermal conductivity while only taking the solid portion into account 


and 


(1 - /3)p v {3p c 

e = e v H e c 

Ps Ps 


__ (1 - P) Pv , (3p c 
cx — 0 l v -r ol c 

Ps Ps 


(34) 


(35) 


II.E. Porous Flow Law 

Applying the Navier-Stokes momentum equations to flow through the char layer would require detailed knowledge 
of the pore structure, and that information is typically not known. Consequently, a porous flow law can be used 
as simplified momentum equations. Porous flow laws typically require extra knowledge about the material beyond 
thermophysical properties. These properties can include the porosity and permeability of the solid, as well as the 
viscosity of the gas flowing through the porous medium. The porosity and permeability can be determined through 
material testing and are provided to the code as a function of extent of reaction. For development purposes, Darcy’s 
law 26 is being used to govern the porous flow. Through experiments, Darcy determined how the volumetric flow rate 
of a laminar flowing fluid is related to the pressure gradient within the fully saturated porous medium. 


Q = -A-VP 
P 


(36) 


The superficial or filtration velocity is the volumetric flow rate averaged over the cross-sectional area of the medium 
and is given by 


=5 = -“ V p 


A 


P 


( 37 ) 
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The average or seepage velocity of the fluid is the volumetric flow rate averaged over the cross-sectional area through 
which the fluid can flow (porous area) and is given by 


v 


9 


Q_ 

<j>A 


-fvp 


(38) 


which assumes that the surface porosity is equal to the volumetric porosity. Darcy’s law is valid for steady laminar 
flows with “sufficiently” low Reynolds numbers. While it is conceded that the pyrolysis gas flow within the porous 
ablator may not be adequately described by Darcy’s law for many situations, it is being used for code development 
purposes due to its simplicity. The addition of more appropriate physical models can certainly be accomplished once 
sufficient data is available to support such a change and the necessary material properties are known. 

The seepage velocity can be used to determine the gas mass flux (mass flow rate per unit area) at any point within 
the medium according to 


ling = (<t>p g ) Vg = - O Pg) ^ V. P = ~Pg^VP = P gV' g (39) 

Each individual velocity component can be calculated according to 

= ( 4 °) 

where k denotes the component and k,/ c denotes the k th row vector of the permeability tensor. Substituting Darcy’s 
law into Eq. (3) gives 

dt ^=u g + V ^p g *VP^+V-(<j>p g v m ) (41) 

Notice that Darcy’s law models the convective term as a diffusive term analogous to Fourier’s law for heat conduction 
and Fick’s law for diffusion. For stationary mesh problems, this results in a character change from hyperbolic to 
parabolic which aids in the stability of the solution procedure. 


ILF. Galerkin Weak Statement 


The governing equations can be restated as 


and 


where 


Energy: - V • (fcVr) + V • (, <j>p g h 0g v g ) - V • (ph Q v m ) -Q = 0 


Q p s 

Solid Mass: — - — u> s — V- (p s v m ) = 0 
at 


Gas Mass: 


dt 


LUg + V ■ (<l>PgVg) - V • (<J bpgV m ) = 0 


Darcy: v g 


-fvp 


(42) 

(43) 

(44) 

(45) 


Since the solid mass equation on a stationary grid only depends on the temperature and density history of a given 
point, it is acceptable to solve the solid mass equation on a stationary mesh while solving the other equations on a 
moving mesh. Consequently, the solution procedure will be to solve for the temperature and gas density fields through 
integration of the partial differential equations (PDEs) given in Eqs. (42) and (44) while the solid density will be treated 
as a nonlinearity, and the resulting ordinary differential equation (ODE) governing the solid density evolution will be 
numerically integrated. 


II.F.l. Energy Equation 

A Galerkin weak statement can be developed for the energy equation (Eq. (42)) by first multiplying it by a suitable 
test function, v, and integrating over the domain, O with boundary, T, while integrating the second, third, and fourth 
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terms by parts to give the natural boundary condition terms. The weak statement is then: Find T £ H 1 such that 


In L 


d{pe 0 ) , ^ 

v-r- h Vu- 

at 


dTl 


(kVT^j - Vt; • (4>p g h 0g Vg) + Vr • ( ph 0 v m ) - vQ 

+ j) (yh 0g m g - vh 0 m m + vq conds ) dT = 0 Vu <E P,) (46) 


where the boundary mass flux due to gas convection is 

rhg = {(j>p g ) v g ■ h 


(47) 


the boundary mass flux due to mesh motion is 


TTlm — P^m ' ^ 


(48) 


II.F.2. Gas Mass Conservation Equation 

Likewise, a Galerkin weak statement can be developed for the gas mass conservation equation.Eq. (44). Again, the 
equation will be multiplied by a suitable test function and integrated over the domain while integrating the third and 
fourth terms by parts to give the natural boundary condition terms. The weak statement is then: Find (j)p g £ // 1 such 
that 

J ^ V ' p 3 v g) + ^ v ■ (fiPgVm) + rn s v'j dfl + j> [u ( m g - fn gm )] dT = 0 \/v £ Hq (49) 

where the pyrolysis gas mass flux at the boundary due to mesh motion is 

% = i&Pg) v m ■ n (50) 


II.G. Finite Element Formulation 


Eqs. (46) and (49) can be discretized by expanding the independent variables and test functions in terms of a finite 
dimensional basis 


N 


N 


( peo) h (®) = X! (P e °h (*) ( pho) h (at) = ^ (pMj i’j ( x ) 

7=1 7=1 

N N 

(' <t>Pg)h (*) = (M) 7 i’o ( x ) T h (at) = Y T 7^7 (*) 

7=1 7=1 

N N 

VT h (at) = £ Tj V ipj (at) VP fc (at) = ^ P, V0,- (at) 

i=i i=i 

A/ - AT 

Ufc (at) = ^ Vilpi (at) Vr k (at) = Y ^ (at) 


i=l 




(51) 


where the subscript h is introduced to denote a finite dimensional approximation. In this study, first order Lagrangian 
basis functions are used to give a second order accurate spatial discretization, and the verification of this order of accu- 
racy will be demonstrated in Section III. Notice that pe Q , ph 0 , and T behave according to the same finite dimensional 
basis. While this may introduce physical inconsistencies within and element, the errors vanish as the mesh is refined. 
This was chosen to simplify implementation and to improve computational efficiency. This assumption will be shown 
to not affect the spatial order of accuracy of the scheme. During the solution process, pe 0 and ph 0 will be evaluated at 
the nodes with the known nodal temperatures and densities, and all four quantities will be interpolated to the element 
interior using the basis functions. An analogous situation exists for the pressure gradient with respect to temperature 
and density. 

Since the unknowns are no longer functions of x , the PDE system reduces to an ODE system in which the temporal 
derivatives can be defined as 
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d(pe 0 ) h 

dt 

d(4 , Pg)h 

dt 


3 = 1 

N d^p g ) 3 


E 


dt 


i’j (*) 


(52) 

(53) 


i= i 

Since the equation system should be satisfied for all combinations of nodal shape function coefficients, 17 , their choice 
is arbitrary as long as a unique combination is chosen for each node . 27 In typical finite element fashion, the coefficients 
are chosen to be v, = da for the I th nodal equation. Consequently, the nodal residual equations resulting from Eqs. (46) 
and (49) can now be written 


f ip i ^^ A dn+ f Vt/’t • kVT h dn - f (<t> Pg ) h h 0g v^i-v g dfi 

In at J n J n 

+ f (ph 0 ) h Vil’i ■ v m dfl + [ tpiq C ond 3 dT + [ ipih 0g m g dT - [ ipih 0 m m dT - [ ipiQdfl = 0 (54) 


In 


in 


and 


, dt 


dfi — {<t>p g ) h V0i • v g dn+ / {<t>p g ) h Vt/’i • v m dn 

Jo Jo 


+ / 'i/ji'rhgdT — / ipiTTig^dr + / 'ip i m s dtt = 0 (55) 


for i = 1,2, ...,N. 


II.H. Time Discretization 


The semidiscrete weak form of the system given by Eqs. (54) and (55) is discretized in time using backwards finite 
difference schemes. Both first and second-order accurate schemes may be derived from Taylor series expansions in 
time about U h (t„+i) = U n+1 : 


U n = U 


n+1 


(<„ - «„+.) + ‘" +l)2 + o ((<„ - („ +1 ) 3 ) 


U n -i — U n+X + 


dU n+1 
dt 

dU n+l . . d 2 U n +\ (t„_i — tn+l ) 2 


(in-1 — in+l) + 


dt ~“ _r±/ 1 dt 2 2 

These expressions can be manipulated to create difference formulas of the form 28 - 29 

dU n+ 1 


+ O 


dt 


— OltU n+l + PtU n + 'JtUn-1 + O (At^ +1 ) 


(56) 


to yield either a first or second-order accurate scheme. The weights a t , dt, and 7 1 are given for p = 1 and p = 2 in 
Table 1. 


Table 1. First and second-order accurate time discretization coefficients. 


p 

OL t 

Pt 

It 

1 

1 

-1 

0 

9 


1 l 1 

Atn + l 


Pt it 

At n _|_ x 1 A t„_ 

At„(Ai n+1 +At n ) 


The resulting temporal derivatives are 

= ot t (pe D )” +1 + f3t (pe 0 )j + 7 1 ( P e o )” _1 (57) 

= a t (PPg)™ +1 + Pt (■ 4>Pg) n j + 7 1 (58) 
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ILL Linearization with the Complex-Step Method 

Eqs. (54) and (55) form a discrete nonlinear system of equations with unknown nodal values T n+1 and (0p g )” +1 . 
This nonlinear system will be solved with a typical Newton scheme. Since the intent is to reach a fully implicit 
formulation, it is necessary to linearize the residual equations in iteration space. This will be done according to the 
familiar Taylor-series expansion 


7er +1 = n 



AT) 


d Ki 


d(<t>Pg)i 


A (4>Pg) 


gjj 


higher order terms 


(59) 


where the superscript u has been introduced to denote iteration level and 


ATj : TJ +1 - Tj 


(60) 


and 


A (<t>Pg)j = ( <t>Pg)] +1 - {<t>PgTj 


(61) 


Keep in mind that the Jacobian terms, jyy, are only nonzero if the support of the basis functions ti> t and yjj overlap. 

Since the intent of the code is to allow for easy addition of material models and porous flow laws, the Jacobians 
will not be derived analytically, rather they will be calculated exactly (to machine precision) using the complex-step 
method. The implementation of the complex-step method does not disallow the use of analytically derived Jacobian 
terms (or some mix of the two), so these terms can be replaced with analytical expressions at a later date if desired. 

The residuals can be expanded according to a Taylor-series expansion in independent variable space about a point 


Tj-, ( < t > Pg)j 


. For example, an expansion with complex step in temperature, iAT, is 


n 


Tj + iATj, (<pp g ) 


gjj 


= n 


T j’(<t>Pg) 


gjj 


+i w, {ATs) ~ 


ia2 V 


2 d Tj 


6 DTf 


higher order terms 

(62) 


Now looking at just the imaginary part and ignoring the higher order terms gives 


SfiK 


Tj +iAT j ,(<j)p g ) j | 




dTj 


6 d T? 


Solving for the first derivative gives the Jacobian terms 


on 




Tj + iATj,(<f>p g )j | 

ATj 


+ o 


(AT# 


(63) 


(64) 


This process can be repeated for each independent variable. The proper choice for the perturbation steps is not 
immediately obvious from looking at the equations. In general, the perturbation step should have some relationship to 
the order of magnitude of the independent variable (in this case Tj). 


ATj = rTj 


(65) 


Taking advantage of the fact that the problem will be solved on a finite precision machine, r can be chosen so that 
the second order error term in the Jacobian expression in Eq. (64) is smaller than the smallest order represented in 
the first term. Consequently, by taking advantage of truncation error, the derivatives will be accurate to machine 
precision. Having such a small step does not affect the division operation in the first term of Eq. (64) because multi- 
plication/division operations simply result in an exponent shift, r = 10" 20 is used for the calculations in this study, 
and an investigation to the sensitivity of this parameter will be presented in Section IV. 

A complex residual calculation is required for each Jacobian term. The scope of the residual recalculations can be 
reduced by only recalculating those terms that depend on the degree of freedom with the current complex perturbation. 
One fortuitous aspect of this method is that there is no separate residual calculation required for the unperturbed state 
(Right Hand Side of linear system). If only one independent variables is perturbed (use T for example), then the real 
part of Eq. (62) can be solved for the residual 


n 


Tj’ ( ( t > Pg)j 




T, 


iATj, (op g )j 


O 


(AT,) 2 


( 66 ) 
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Again since the independent variable perturbation has been chosen to be sufficiently small, the order of magnitude 
of the error terms is smaller than the smallest order represented in the complex residual calculation. This results in 
a residual calculation exact to machine precision. Through calculation of the Jacobian terms with the complex step 
method, the calculation of the real residual comes free. 

The advantage of the complex-step method is in the development cost. The derivation, implementation, and 
debugging of analytical Jacobian terms can be time consuming. Plus, the addition and/or modification of models 
with the complex-step method can be accomplished with relative ease. Additionally, most codes with complicated 
governing equations employ approximate Jacobians which will result in less the optimal convergence of the nonlinear 
solution routine. The drawback of the complex-step method is the extra time spent doing complex arithmetic. While 
more memory is required to store complex variables, the necessity for complex variables only occurs at the element 
assembly level. Therefore the memory requirements for the complex-step method are insignificant relative to the 
memory requirements for storing the mesh, variables at several time levels, and the linear system. 

II J. Additional Numerical Issues 

The integrals in Eqs. (54) and (55) are evaluated numerically with using Gaussian quadrature, and the resulting linear 
system is solved using the Generalized Minimal Residual method with preconditioning. 21 The code was developed 
using the libMesh finite element library. 30 Essentially, the library handles all of the physics independent aspects of 
the code including I/O, mesh connectivity and element data, domain decomposition, parallel communication, adaptive 
mesh refinement, and linear system solution. libMesh provides a seamless interface to several other libraries including 
METIS 31 for domain decomposition and PETSc 32 for solving the linear system. 

ILK. Boundary Conditions 

The code can handle a rich suite of boundary conditions, including multiple flux type boundary conditions that can be 
summed over a given boundary. Only those that are used in this study will be described here. 


II.K.l. Specified Temperature 

Dirichlet conditions and other conditions that require the substitution of a nodal residual equation with a boundary 
condition specific equation will be imposed via the penalty boundary method (PBM). 33 With this method, the L 2 
projection of the residuals are added to the matrix. The projection is multiplied by some large factor so that, in floating 
point arithmetic, the existing (smaller) entries in the matrix and right-hand-side are effectively ignored. The boundary 
integral in the weak form becomes 

- [ ip t lldr = 0 (67) 

e Jr 

where eCl. For the specified temperature condition, the residual equation is given by 

Tt = T spec - T h (x) = 0 (68) 

where T spec is the known specified temperature. Substituting in the definition of Th (x), Eq. (51), gives 

N 

TZ = T spec - ^2 ( x ) = 0 (69) 

3 = i 

Substituting the residual definition into Eq. (67) gives the weak form of the discrete residual equations 


TZ, = ^ J ipi {^Tspec - X] T /Vb dT = 0 


(70) 


Employing a numerical integration technique, such as Gaussian quadrature, to evaluate the surface integral gives 


# QPs 


*< = 7 £ 


QP = 1 


N 


J xWipi ( T spec - ^2 Tjipj 
j = i 


(71) 


QP 
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where QP denotes the quadrature point and J x VI' is the Jacobian of the Uansformation matrix times the integration 
weight. Linearizing in iteration space according to a Taylor-series expansion while ignoring higher order terms gives 


n 


v+l _ V v , VP-t 

1 + 2^ frr, 


j = i 


A Ti 


Differentiating Eq. (71) gives the Jacobian terms for the linear system 

dPi JxW 


dTi 


-tpiipj 


(72) 


(73) 


Given the ease of the Jacobian derivation, this boundary condition has been implemented using analytical Jacobians 
as opposed to using the complex perturbation approach. 


II.K.2. Specified Heat Flux 

The specified flux boundary condition is linear and requires no Jacobian contributions. The boundary heat flux term 
in Eq. (54) becomes 

l IpiQcondgdf-' = I tpiq spec dT (74) 

where q spec is the known specified flux. Note that the sign change is intended to simplify usability. In the original 
PDE derivation, fluxes into the body were considered a negative quantity. So the sign change allows the user to 
input a positive heat flux when the intent is to put heat into the body. Introducing Gaussian quadrature for numerical 
integration gives 

„ # QPs 

I tpiq S p ec dP — ^ ' J x W qspec^Pi (75) 

"' r QP = 1 


II.K. 3. Convection 

The convective heat flux is given by 

Qconv 77 (T r Tyfi) 

Substituting into the boundary heat flux term in Eq. (54) gives 


J ipiQ con d s dT = - J tpiH (T r - T w ) dT 


(76) 


(77) 


Introducing Gaussian quadrature for numerical integration and the finite dimensional basis representation of tempera- 
ture gives 


, #QPs 

- / i\H (T r - T UI ) dT = 

QP= 1 


N 


-J x WtfiH ( T r - 

j = i 


(78) 


QP 


The Jacobian contributions for the convective boundary condition are calculated using the complex perturbation 
method. 


II.KA. Specified Gas Density 

The specified gas density condition is a Dirichlet condition on the gas continuity equation. In real world problems, the 
surface density is not likely known. Therefore it is expected that this condition will be infrequently used. However, 
it is very useful for debugging purposes and its implementation is verified in this study. The specified gas density 
condition is imposed via PBM, and the residual equation is given by 

K = 0 i>P 9 )s P e C (0 P 9 )h C * ) = 0 (79) 

where ( 4>Pg) spec is the known specified density. Substituting in the definition of (( t>P g ) h { x )> Efi- (51), gives 

N 

n = (M) S pec - ( < M?)j V’J (*) = 0 (80) 

3 = 1 
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Substituting the residual definition into Eq. (67) gives the weak form of the discrete residual equations 


7 Zi = 




N \ 

Pa)j i’j O) Ur = o 


( 81 ) 


Employing a numerical integration technique, such as Gaussian quadrature, to evaluate the surface integral gives 


# QPs 


** = 7 £ 


QP—1 


N 


J x Wipi (# ff ) spec -E {<PPg)j i’j 


j = 1 


I QP 


(82) 


and differentiating Eq. (82) gives the Jacobian terms for the linear system 


dfcj 

d(<t>Pg)j 


JxW 

e 




(83) 


Given the ease of the Jacobian derivation, this boundary condition has been implemented using analytical Jacobians as 
opposed to using the complex perturbation approach. Although the gas density degree of freedom is specified with this 
condition, the energy equation still has a pyrolysis gas mass flux term that must be evaluated along with its Jacobians. 


II.K. 5. Specified Mass Flux 

The specified mass flux condition effects terms in both the energy and gas continuity equations. 


and 


/ ^ ih 0g rn spec dr 
r 

(84) 

J fiirhspecdT 

(85) 


Like the specified heat flux condition, the negative sign is introduced so that the user inputs a positive mass flux to 
denote flow into the body. This condition is linear for the gas continuity equation, however it is nonlinear for the 
energy equation since the gas total enthalpy depends on both density and temperature. Consequently, Jacobian terms 
are required, and these again are calculated through the complex-perturbation method. 


II.K. 6. Specified Pressure 

Since the independent variables in the linear system are temperature and gas density, the specified pressure condition 
is nonlinear and must be applied through PBM to replace the appropriate gas density residual equations in the linear 
system. The residual equation for the specified pressure condition is 

P = Pspec - P(x) = 0 ( 86 ) 

where P spec is the known specified pressure. Substituting the residual into Eq. (67) gives the weak form of the discrete 
residual equations 

7 li = -[ ipi ( Pspec - P) dT = 0 (87) 

e Jr 

Again using Gaussian quadrature to evaluate the surface integral gives 

# QPs 

7 li = - J2 {JxWA (Pspec -P)]\ QP (88) 

6 QP—l 

The Jacobian entries are calculated using the complex-perturbation method. The energy equation still has a pyrolysis 
gas mass flux term that must be evaluated along with its Jacobians. 
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III. Verification 


Verification exercises are typically done by generating solutions on successively refined space or time meshes for 
problems with known exact solutions. Examining the solution error as the mesh is refined shows the order of accuracy 
of the model as implemented. Typically, the spatial and temporal orders of accuracy are independently verified. If the 
discretization error is of the form 


DE = ciAt p + bAh q + (higher order terms) 


(89) 


Then the temporal order of accuracy can be determined by making bAh q negligible through the use of a very fine mesh 
and solving a problem with successively refined time steps. The logarithm of the discretization error as a function of 
the logarithm of the time step should have a slope of p 

g[log (DE)] = 

<9[log (At)] P 

with intercept log (a). Similarly for the spatial order of accuracy, aAt p can be made negligibly small through the use 
of an extremely small time step. The logarithm of the discretization error as a function of the logarithm of the mesh 
spacing should have a slope of q 

<9[log (DE)] _ 

9 [log (Ah)] q 

with intercept log (b). This method can be time consuming considering that the problems either need to be solved with 
an extremely fine mesh or time step. 

In the current study, an alternative approach is considered in which the spatial and temporal orders of accuracy are 
simultaneously verified. 15 Factoring out a bAh q from the discretization error gives 

/ nAt Pth \ 

DE = 1 + — bAh qth (90) 

V bAh^h J 


where the subscript th is introduced to denote the theoretical orders of accuracy. If it is assumed that the theoretical 
orders of accuracy p t h and q t h are known, which they usually are, then simultaneous refinement of the spatial and 
temporal meshes while holding bAh^th constant can be performed. If the observed slope of the log-log relationship 
between the discretization error and mesh spacing, q Q b s , is the same as the expected theoretical slope, q t h. 


0[log (DE)] _ 

3 [log (Ah)] ~ Qth ~ q ° bs 


(91) 


with intercept 


log 


b 


( a A t Pth \ 

\ bAh qth J 


then both the temporal and spatial orders of accuracy are simultaneously confirmed as long as ' s on the order of 
or greater than 1 . This is because one would like the intercept to be sensitive to this quantity so that a non-co-linear 
result would suggest that either p Q b s or q„/, s is not what is expected from theory. The drawback of this methodology 
is that if the expected result is not produced, then the developer cannot be certain whether the error is in the temporal 
or spatial discretization. However, if the implementation of the temporal discretization has been previously verified, 
then the error would have to lie in the spatial discretization. This method offers a significant time savings over the 
traditional independent refinement approach since it is not necessary to solve the problem with extremely small time 
steps or fine meshes. 

When comparing a numerical solution to an exact solution there are several different acceptable forms of the 
discretization error metric. These could include the value of independent variables at a specific location and/or time, 
or the metric could be global comparison of the solutions. For this study, the discrete normalized root-mean squared 
(RMS) error of a nodal quantity of interest across the domain is chosen as the error metric. 


DErms 


1 

Qo\ 


i N 

n—1 


(Q exact Q sim) 


( 92 ) 
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where N denotes the number of nodes, exact denotes the exact solution, sim denotes the simulated solution, Q denotes 
the quantity of interest (usually the independent variable(s)), and Q a represents the normalization value (usually the 
initial value). When evaluating the error as a function of mesh refinement, the grid metric has been chosen to be 

Number of Elements in One Dimension: d ' r \[W e 


where dim is the number of spatial dimensions in the problem and N e is the total number of elements in the domain. 
The number of elements in one dimension is inversely proportional to the characteristic mesh length scale. Ah. Since 
for a given problem, the characteristic domain length, L, is constant, the characteristic mesh length scale can be defined 
as 



(93) 


Substituting this relationship into Eq. (90) gives 


DErms 


This leads to the verification criteria 


aAt Pth 

bAh qth 




9 [log (DErms)] 

d[log ( d W)] ~~ 9th ~ ^ 


(94) 


(95) 


This means that in order to simultaneously confirm the spatial and temporal orders of accuracy, a slope of — q th in 
log {DErms)- log ( di \fN^) space must be observed when solving the problems on successively refined space and 
time meshes while holding constant. 

In this study, second order spatial accuracy ( q = 2) and second order temporal accuracy (p = 2) will be verified 
through the simultaneous space/time refinement approach outlined above. Additionally, it is the intent to verify both 
the two and three dimensional capabilities for all allowable element types. The element types included in this study 
are: triangles and quadrilaterals in two dimensions and hexahedrons, tetrahedrons, and prisms in three dimensions. 
While verification studies have been performed with pyramids, none of the problems included in this study were solved 
with pyramid elements. 

In addition to verifying the temporal and spatial orders of accuracy, it is necessary to check the rate of convergence 
for the Newton scheme. Theoretically, a correct implementation of Newton’s method should exhibit second-order 
convergence. If the solution converges at some other rate, it is likely that a mistake was made in the derivation 
and/or implementation of the method, likely in the Jacobian terms. The theoretical relationship between error, E, at 
successive iterations is given by 

E v+1 = a(E v ) 2 (96) 

In this study, the relative L 2 of the independent variable correction vector is used as the error metric. 


E v 


(l 2 Y 

( L 2 )" =1 


(97) 


When plotting the logarithm of the error at a given iterate, v + 1, versus the logarithm of the error at the previous 
iterate, v , the function should be linear with a slope of 2 and intercept log (a). 


9 [log E v+1 ] 
3 [log E v ] 


(98) 


It is evident that no point will exist for the final (converged) iteration in Newton’s method since the error metric 
is zero. It is also important to note that demonstrating second-order nonlinear convergence does not imply that the 
method is converging to the correct answer for a given problem. Similarly, the absence of second-order convergence 
does not imply that the solution is incorrect. Second-order nonlinear convergence simply proves that the Jacobian 
matrix contains the correct values for the implemented residual equations, and the efficiency of the solution procedure 
is maximized. Incorrect residual equations can still exhibit second-order convergence if the associated sensitivities 
were correctly derived and implemented. 

The verification problems outlined in this section are a small subset of a large suite of problems. They present 
the general strategy for verification that has been adopted in this code development process. In general, a verification 
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problem exists for each boundary condition, and each exercises several terms in the governing equations. Currently, 
there is no verification problem that exercises all nonlinear terms in the coupled governing equations at the same time, 
but this is the subject of future work. The verification problems also serve as the regression test suite. Verification 
and regression tests are performed across multiple architectures including compilers and MPI libraries. Finally, all 
solutions presented in this study were solved in parallel, but the verification and regression tests are also performed in 
serial. 


Ill A. Transient Conduction with Variable Properties 

III.A.l. Problem Statement 

This problem has been previously used to verify a thermal response code as outlined by Amar, 15 and it is intended to 
verify the implementation of the two dimensional nonlinear transient heat conduction terms in the governing equations 
and their requisite terms for Newton’s method. Consider a non-ablating, non-decomposing, two dimensional, uniform 
density planar slab of equal length and width subject to a constant specified heat flux on one face, x = 0, and adiabatic 
on the other three faces. Since conduction will only be occurring in one dimension, the governing equation, boundary 
conditions, and initial conditions that characterize the problem are 


pC p 


dT 

~9t 


dx \ dx ) 


for 0 < x < L 



— Qspec 

x—0 


-k 


dT 

dx 


= 0 

x—L 


> 


(99) 


T(x,t = 0) = T a 


and the parameters for this problem are given in Table 2. The thermal conductivity and specific heat are given in 
Table 3 and are assumed to vary linearly with temperature while the thermal diffusivity is constant giving 

k(T) = k 1 +*^^r(T-T 1 ) (100) 

C p (T) = C Pt i + Cp T 2 ~ ° pA (T - Ti) 

±2 — ±1 

__ hi _ fc 2 _ k ( T ) 

pC p ,i pC Pt 2 pCp (' T ) 


( 101 ) 

( 102 ) 


Table 2. Variable properties problem parameters. 


T 0 = 300 K 


g spec = 7.5 x 10 5 W/ m 2 

L = 0.01m 

p = 8000 kg/ m 3 

= k/{pC p ) = 2.5 x 10 -6 “ 2 /s 


Table 3. Material properties for temperature dependent properties verification problem. 


T,K 

k, w/ m . K 

Cp, J /kg K 

= 300 

10 

500 

T 2 = 1300 

100 

5000 
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III.A.2. Analytic Solution 

The analytical solution to this problem can be found by first applying Kirchhoff’s transformation as described by 
Arpaci 34 and then using the related constant property solution shown in Beck et al. 35 . The transformed energy equation 
now has the same form as the related constant property case, and the analytical solution is 


where 


6 (x, t) — 9 0 at 1 

QspecL / L“ 3 



9 = 


(T-T 1 ) + 


&2 — k\ 1 
T 2 - Ti 2k x 


(T-Ti) 2 


(103) 


(104) 


III.A.3. Grid Refinement and Nonlinear Convergence Studies 

Although the solution is governed by a one dimensional differential equation, the problem can be used to verify the 
two dimensional heat conduction capability in the code. Typically, when verifying two or three dimensional codes, the 
developer would like to exercise flux terms in all directions. This is especially important when dealing with solvers 
that only use structured meshes since fluxes in the i , j, and k directions may be evaluated with different routines. For 
an unstructured code, fluxes in all direction are not necessary. This is especially evident when dealing with triangular 
elements since no faces on the interior of the domain are aligned with the flux direction. The problem was solved on a 
series of four quadrilateral and four triangle meshes as outlined in Tables 4 and 5 respectively, and the coarse meshes 
can be seen in Figure 2. Each mesh was uniformly refined by a factor of two to reach the next grid level. 

Table 4. Quadrilateral mesh parameters for variable property conduction problem. 


Grid Name 

N 

N e 

VNe 

Ai, sec 

Coarse 

441 

400 

20 

4 

Medium 

1681 

1600 

40 

2 

Fine 

6561 

6400 

80 

1 

Very Fine 

25,921 

25,600 

160 

0.5 


Table 5. Triangle mesh parameters for variable property conduction problem. 


Grid Name 

N 

N e 

VNe 

At, sec 

Coarse 

787 

1492 

38.6 

4 

Medium 

2980 

5798 

76.1 

2 

Fine 

11,757 

23,192 

152 

1 

Very Fine 

44,322 

88,002 

297 

0.5 


The grid convergence results can be seen in Figure 3. For this problem, the normalized L 2 of the temperature 
solution error vector is used as the error metric where Q 0 = T a = 300 K in Eq. (92). It is evident that both element 
types exhibit second order spatial and temporal accuracy. Also as expected, the quadrilateral elements provide a 
significant accuracy advantage over triangles. 

Since this problem is nonlinear, it is also necessary to verify the convergence rate of the Newton scheme. Fig- 
ure 4 shows the nonlinear convergence results for the first time step on the coarse quadrilateral mesh. This result is 
representative of the convergence history throughout the problem for any of the meshes in this study. 

III.B. Three Dimensional Transient Conduction with Perfect Gas Convective Boundary Conditions 

III.B.l. Problem Statement 

Consider a Cartesian aligned uniform density constant property cube whose center is at the origin and has side length, 
2 L, that is exposed to a convective heat flux with constant recovery temperature on all faces, but the heat transfer 
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Quadrilateral 


Triangle 


Figure 2. Coarse two dimensional meshes. 



Number of Elements in One Dimension 


Figure 3. Grid convergence results: variable property conduction problem. 
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Figure 4. Nonlinear convergence results: variable property conduction problem. 
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coefficients are only the same on opposing faces of the cube. The governing equation, boundary conditions, and initial 
conditions that characterize the problem are 


pC'p— — V-fcVT = 0 

fcVT • h = H x ( T r — T) for x = —L and x = L faces 

fcVT • h = H y (T r — T) for y = —L and y = L faces * 

k'VT ■ h = H z (T r — T) for z = —L and z = L faces 

T(x,t = 0) = T a 

and the parameters for this problem are given in Table 6. 

Table 6. Cube convection problem parameters. 


(105) 


T 0 = 300 K 

T r = 1300 K 
L = 0.005 m 

p - 8000 k s/m 3 

C p = 500 J/kg-K 
k = 10 w/ m . K 

a = k/(pC p ) = 2.5 x 10~ 6 m 2 / s 
H x = 1500 w/ m 2 K 
H v = 1000 w/ m 2 K 
H z = 500 w/ m 2 K 


III.B.2. Analytic Solution 

One dimensional solutions to the heat equation with constant thermal conductivity and specific heat can be combined 
to give exact solutions to simple two and three dimensional geometries as outlined in Bejan. 36 For this problem, the 
exact solution used is for the one dimensional constant property planar slab of length L with front face convective 
environment and an adiabatic back face as outlined in Arpaci. 34 The exact solution in the a:-di recti on is given by 


Ox (x, t) 


T( x, t) — T r 
T — T 


2E 

n — 1 


sm v n 


v n + sin v n cos v n 


exp (—v^at/L 2 ) cos (v n x/L) 


(106) 


for all u n > 0 where the eigenvalues can be found by solving the implicit relationship 


v n sin u„ = Bi cos v n 


(107) 


and the Biot number is given by 



(108) 


and there is an analogous expression for 9 y (y. t) and 9 Z (z, t ). The one dimensional solutions can be combined to give 
the three dimensional distribution according to 


T(x,t.) - T r 
T 0 — T r 


9 x 9 y 9 z 


(109) 


The above expression is only valid for the octant where x > 0, y > 0, and z > 0. Since the solution is symmetric 
across the x = 0, y = 0, and 2 = 0 planes, Eq. (109) can be appropriately applied to give the temperature distribution 
in the entire domain. The solution after 20 seconds on the cube’s exposed faces can be seen in Figure 5. 


22 of 38 


American Institute of Aeronautics and Astronautics 





Figure 5. Solution at 20 seconds for cube convection problem. 


III.B.3. Grid Refinement Study 

The problem was solved on a series of three hexahedron, three prism, and three tetrahedron meshes as outlined in 
Tables 7, 8, and 9 respectively, and the coarse meshes can be seen in Figure 6. For this problem, the normalized L 2 of 
the temperature solution error vector is used as the error metric where Q a = T a = 300 K in Eq. (92). The results of 
the grid refinement study can be seen in Figure 7, and it is evident that the code exhibits both second order spatial and 
temporal accuracy. As expected, the hexahedrons were the most accurate for a given mesh size while the tetrahedrons 
were the least accurate. Although the boundary condition is nonlinear, the interior solution is linear resulting in a very 
weakly nonlinear problem. Consequently, the Newton scheme converged in too few iterations to generate a meaningful 
trend for a nonlinear convergence study. 

Table 7. Hexahedron mesh parameters for cube convection problem. 


Grid Name 

N 

N e 


At, sec 

Coarse 

9261 

8000 

20 

0.5 

Medium 

68,921 

64,000 

40 

0.25 

Fine 

531,441 

512,000 

80 

0.125 


Table 8. Prism mesh parameters for cube convection problem. 


Grid Name 

N 

N e 


At, sec 

Coarse 

16,275 

29,360 

30.8 

0.5 

Medium 

122,754 

233,040 

61.5 

0.25 

Fine 

956,286 

1,863,200 

123 

0.125 
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Table 9. Tetrahedron mesh parameters for cube convection problem. 


Grid Name 

N 

N e 

m 

At, sec 

Coarse 

19,414 

103,268 

46.9 

0.5 

Medium 

130,932 

736,834 

90.3 

0.25 

Fine 

957,152 

5,571,020 

177 

0.125 



Figure 6. Coarse cube meshes. 
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Figure 7. Grid convergence results: cube convection problem. 
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III.C. Three Dimensional Transient Conduction with Manufactured Solution 

Often it is difficult to find an analytic solution that exercises more than a few terms in the governing equations of 
interest. It is desirable to have an exact solution to a problem that exercises not only the nonlinear terms in the govern- 
ing equations but also several code options simultaneously. One common approach for developing such a verification 
problem is the Method of Manufactured Solutions. 20 22 For a general PDE with temperature as the independent vari- 
able 

V{T) = 0 (110) 

where V represents the differential operator, a simplified outline to the approach is as follows: 

1. Choose the desired solution, T ( x , t), to the PDE. 

• Solution should be composed of simple analytic functions that are functions of both space and time if 
desired. 

• Solution should be sufficiently smooth so that the theoretical order of accuracy can be achieved. 

• Solution should exercise all desired terms in the PDE. 

• Solution should be differentiable at least as many times that is dictated by the PDE. 

• Solution should not violate physics such that the software will not run. For example, negative densities or 
temperatures on an absolute scale should not be part of the solution. 

• Like all verification problems, the solution should not be limited to something that is physically meaningful 
or of practical importance since verification is simply a mathematical exercise. 

2. Choose analytic model equations for any properties that depend on the independent variable (e.g. thermal 
conductivity and specific heat for the energy equation). 

3. Substitute the solution and property models into the PDE, and allow the differential operator, V, to operate on 
these functions. 

4. In general, the chosen solution is not a solution to the PDE, V ( T ) / 0. However, the PDE can be redefined to 
have a source term. 

V(T) = Q(x,t) (111) 

where Q (x,t) is the source term which is defined by the result of step 3. 

5. Implement the resulting source term function in the software since T is the exact solution to the modified PDE 
in Eq. (Ill) with known source term. 

6. Perform verification studies to check for orders of accuracy and convergence rate 


III.C. 1. Solution Development 

The transient temperature distribution, T (x,t), in a constant density material is governed by 

pCp-^r - V-(fcVT) = Q (112) 

and the manufactured temperature solution is chosen to be a function of cosines 

T(x, t) = 9 cos(H x a: + A t t) cos (B y y + B t t ) cos (C z z + C t t) cos (D t t) + T a (113) 

where T a is the average temperature about which the solution oscillates with amplitude 9. The models for thermal 
conductivity and specific heat are chosen to be quadratic in temperature. 


k ( T ) = fc 0 + k\ 
C'p (T) = C p 0 + C P1 


T — T n 


T~T 0 

9 


T — T n 


C.„ 


T — T n 


(114) 

(115) 
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Substituting Eqs. (113) (114) and (1 15) into Eq. (1 12) and solving for the source term gives 

Q (x, t) =p [ C Po + C Pl cos(A x x + A t t) cos (B v y + B t t ) cos (C z z + C t t) cos (D t t) 

+ C P2 cos 2 (A x x + A t t) cos 2 {B v y + B t t) cos 2 {C z z + C t t ) cos 2 {D t t) } 

{ — 9[A t sin( J 4 x x + A t t) co s(B y y + B t t) cos (C z z + C t t) cos (D t t) 

+ B t cos(^4 x a; + A t t) sin (B y y + B t t ) cos (C z z + C t t) cos (D t f) 

+ C t cos (A x x + A t t) cos (B y y + B t t) sin (C z z + C t t) cos (D t t) 

+ D t cos (Ac a: + A t t) cos {B v y + B t t) cos (C z z + C t t) sin (D t t) } } 

— { [k 0 + ki cos{A x x + A t t ) cos {B y y + B t t) cos (C z z + C t t) cos (D t t) 

+ ki cos 2 {A x x + A t t) cos 2 {B v y + B t t) cos 2 (C z z + C t t) cos 2 {D t t) ] 

[—0 ( A 2 + B 2 + C 2 ) cos (A x x + A t t ) cos (B y y + B t t) cos (C z z + C t t) cos(D t <)] 

(116) 

— 6 [ A x sin(A s a: + A t t) co s{B y y + B t t ) cos (C z z + C t t) cos (D t t) 

( — k\A x sin(A x a; + A t t) co s(B y y + B t t ) cos (C z z + C t t) cos (D t t) 

+ 2k2A x sin(^4a;a; + A t t) cos^^rr + A t t) cos 2 {B v y + B t t) cos 2 {C z z + C t t) cos 2 {D t i) ) 

+ By cos^^rr + A t t ) sin (B y y + B t t ) cos (C z z + C t t) cos (D t t) 

( — k±By cos(j4^x + A t t) sin (B v y + B t t) cos (C z z + C t t) cos (D t t) 

+ 2/c 2 i?y cos 2 {A x x + A t t) sin {B y y + B t t) cos (B y y + B t t ) cos 2 (C z z + C t t) cos 2 {D t t) ) 

+ C z cos^^a: + A t t) cos (B v y + B t t ) sin(C z z + C t t ) cos (D t t) 

( — k\C z cos {A x x + A t t) cos {B v y + B t t) sin(C z 2 ; + C t t) cos (D t t) 

+ 2/c 2 C z cos 2 {A x x + A t t) cos 2 {B v y + B t t) s\n{C z z + C t t) cos (C z z + C t t) cos 2 {D t t) ) } 

Notice that the source term turns out to be a complicated function which is typical of MMS. The source terms can 
be derived by hand, but a symbolic manipulator such as Maple 37 or Mathematica 38 can be very helpful with the 
derivation. Note that this manufactured solution can be used for several verification exercises including 

• One Dimensional Conduction (B y = C z = 0) 

• Two Dimensional Conduction (C z = 0) 

• Steady Conduction (. A t = B t = Ct = D t = 0) 

• Linear (C Pl = C P2 = k\ = fc 2 = 0) 

The problem parameters used in this study are given in Table 10. 

Both specified temperature and specified heat flux boundary conditions can be exercised with this problem. For 
the specified temperature condition, the time and quadrature point coordinates are simply plugged into Eq. (113) to 
give the boundary condition. The boundary heat flux at any point can be determined according to 

Qcond, s = —k"VT ■ n = — [ko + k\ cos^^rr + A t t) cos {B v y + B t t) cos (C z z + C t t) cos (D t t) 

+/c 2 cos 2 (A x x + A t t) cos 2 (B y y + B t t) cos 2 (C z z + C t t) cos 2 (D(f)] 

{—9 [^4sin(^4 x a; + A t t) cos {B v y + B t t) cos (C z z + C t t) cos (D t t)n x (117) 

+ B cos (A x x + A t t) sm(B v y + B t t) cos (C z z + C t t) co s(D t t)n y 
+C cos(^4 x x + A t t) cos {B y y + B t t ) sin (C z z + C t t) cos(D t t)n z }} 

where 

h = n x i + riyj + n z k (118) 

1I1.C.2. Grid Refinement and Nonlinear Convergence Studies 

An arbitrary geometry was chosen for this problem such that not all the faces are Cartesian aligned nor are they all flat. 
The problem was solved on a series of three hexahedron meshes and three tetrahedron meshes as outlined in Tables 1 1 
and 12 respectively, and the coarse meshes can be seen in Figure 8. The initial condition and solution to this problem 

27 of 38 


American Institute of Aeronautics and Astronautics 



Table 10. MMS energy problem parameters. 


T 0 = 300 K 

9 = 60 K 

p = 800 k s/m 3 

C P0 = 500J/kg-K 

C P1 = 50 J/kgK 

C P2 = 5 J/kgK 

k 0 = 30 w / mK 

jfei = 3 w/ m . K 

k 2 = 0.3 W /m K 

A x = By = C z = 1 m -1 

A t = B t = C t = 0.1 sec -1 

D t = 0.05 sec -1 


Table 11. Hexahedron mesh parameters for MMS problems. 


Grid Name 

N 

N e 

m 

At, sec 

Coarse 

18.081 

16,000 

25.1 

1.0 

Medium 

136,161 

128,000 

50.4 

0.5 

Fine 

1,056,321 

1,024,000 

101 

0.25 


Table 12. Tetrahedron mesh parameters for MMS problems. 


Grid Name 

N 

N e 

VNe 

At, sec 

Coarse 

20,657 

107,842 

47.6 

1.0 

Medium 

133,029 

736,709 

90.3 

0.5 

Fine 

933,416 

5,375,076 

175 

0.25 
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Tetrahedron 


Hexahedron 


Figure 8. Coarse grids for MMS problems. 




Figure 9. Solution for the MMS energy problem. 


29 of 38 


American Institute of Aeronautics and Astronautics 


can be seen in Figure 9. The problem was solved with a specified flux on the four long sides while the two smaller 
sides had a specified temperature boundary condition. 

For this problem, the normalized L 2 of the temperature solution error vector is used as the error metric where 
Q 0 = T 0 = 300 K in Eq. (92). The results of the verification study can be seen in Figure 10, and it is evident that 
the code exhibits both second order spatial and temporal accuracy. The nonlinear convergence study confirms that the 
Newton scheme is exhibiting second order convergence as seen in Figure 1 1 . 



Number of Elements in One Dimension 


Figure 10. Grid convergence results: MMS energy problem. 


III.D. Three Dimensional Transient Porous Gas Flow with Manufactured Solution 

1II.D. 1. Solution Development 

The transient gas density distribution in a porous medium with Darcian flow is governed by 

(119) 

If it is further assumed that the medium is isothermal, it has uniform isotropic permeability and constant porosity, and 
the gas has constant viscosity and behaves according to the perfect gas law, then the governing equation simplifies to 

*%-^-( P g VP ')=* 9 (120) 

The manufactured density solution is chosen to be a function of cosines 

p g (x, t) = 6cos(A x x + A t t ) cos (B v y + B t t ) cos (C z z + C t t ) cos (D t t) + p 3o (121) 
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Figure 11. Nonlinear convergence results: MMS energy problem. 
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where p 9o is the average density about which the solution oscillates with amplitude 8. The source term is found by 
substituting the manufactured solution, Eq. (121) , into the governing equation, Eq. (120), which gives 

Ug (at, t) = —(j)9 [At sin ( A x x + A t t) cos ( B v y + B t t) cos ( C z z + C t t ) cos ( D t t ) 

+B t cos {A x x + A t t) sin ( B v y + B t t) cos ( C z z + C t t ) cos (D t t) 

+C t cos ( A x x + A t t) cos (B v y + B t t ) sin ( C z z + C t t ) cos (D t t) 

+D t cos ( A x x + A t t ) cos (B y y + B t t ) cos ( C z z + C t t) sin (D t t) } 

kRT 

{ - Pgo d( y A l + By + Cl ) cos ( A x x + A t t) cos (B y y + B t t) cos (C z z + C t t) cos (D t t) (122) 

P 

—0 2 (A 2 + B 2 + C 2 ) cos 2 ( A x x + A t t) cos 2 ( B y y + B t t) cos 2 ( C z z + C t t ) cos 2 (D t t) 

+9 2 cos 2 (D t t) [A 2 sin 2 (A x x + A t t ) cos 2 (B v y + B t t) cos 2 ( C z z + C t t) 

+B 2 cos 2 (A x x + A t t) sin 2 (B v y + B t t ) cos 2 ( C z z + C t t) 

+C 2 cos 2 (. A x x + A t t ) cos 2 ( B v y + B t t) sin 2 ( C z z + C t t ) ] } 

Like the MMS energy problem, this solution can also be used for several verification exercises including 

• One Dimensional Flow ( B y = C z = 0) 

• Two Dimensional Flow (C z = 0) 

• Steady Flow ( A t = B t = C t = D t = 0) 

Unlike the MMS energy problem, this problem can not be made linear due to the pg'V p g term. The parameters used 
in this study can be seen in Table 13. 

Table 13. MMS gas flow problem parameters. 


T = 300 K 

R = 300 J /k g K 
9 = 0.2 kg/ m 3 
p go = 1 kg/m 3 
k=1x 10 -15 m 2 

4 > = 0.2 

/r = 1 x 10 -6 Pa • s 
A x = By = C z = 1 m^ 1 

A t = B t = C t = 0.1 sec -1 
D t = 0.05 sec -1 


Specified density, specified pressure, and specified mass flux boundary conditions can be exercised with this prob- 
lem. For the specified density condition, the time and quadrature point coordinates are simply plugged into Eq. (121) 
and multiplied by the porosity to give the boundary condition. The boundary pressure at any point can be determined 
by first calculating the density according to Eq. (121) then using the perfect gas law. The mass flux at the boundary 
can be calculated according to 




<j>PgV g h = [ 9cos{A x x + A t t) co s(B y y + B t t ) cos {C z z + C t t ) cos (D t t) + p 9o ] 

P 

{—9 [A x sin^^a: + A t t) cos {B v y + B t t) cos (C z z + C t t) cos (D t t)n x 
+ By cos (A x x + A t t) sin (B v y + B t t ) cos (C z z + C t t) cos (D t t)n y 
+ C z cos(^4a;a; + A t t) cos {B y y + B t t ) sin (C z z + C t t) cos (D t t)n y ] } 


(123) 


III.D.2. Grid Refinement and Nonlinear Convergence Studies 

The geometry used for this problems is the same that was used in the MMS energy problem in Section III.C. The 
hexahedron and tetrahedron grid parameters can be found in Tables 11 and 12 respectively, and the coarse meshes can 
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be seen in Figure 8. The initial condition and solution to this problem can be seen in Figure 12. The problem was 
solved with mass flux, density, and pressure boundary conditions, each on two sides. 




0.99 

0.98 

0.97 

0.96 

0.95 

0.94 

0.93 

0.92 


Figure 12. Solution for the MMS gas flow problem. 


For this problem, the normalized L 2 of the gas density solution error vector is used as the error metric where 
Q 0 = p 9o = 1.0 k s/m 3 in Eq. (92). The results of the verification study can be seen in Figure 13, and it is evident that 
the code exhibits both second order spatial and temporal accuracy. The nonlinear convergence study shows that the 
Newton scheme is exhibiting second order convergence as seen in Figure 14. 


IV. Complex-Step Size Study 

In Section II. I, the relative perturbation factor, r in Eq. (65), for the complex-step method was discussed. In order 
to have a Jacobian exact to machine precision, r should be chosen such that underflow occurs on the right-hand-side 
of Eq. (64). However, this is impossible if the order of magnitude of the Jacobian is not known before hand. In this 
code development project, r = 10 -20 is used. If the independent variable is on the order of 10 3 , which is typical 
of temperatures in a re-entry ablation problem, then Eq. (65) says that the step size will be on the order of 10 -1 ', 
which gives an error term on the order of 1 0 31 in Eq. (64). Consequently, for a calculation with 16 digits of precision 
(typical for double precision), underflow will occur for Jacobians as small 10 -18 resulting in exact Jacobians. Once 
the Jacobian gets smaller it will be approximate, but Jacobians less than 10“ 18 are typically inconsequential to the 
solution of the linear system. Although r = 10~ 2Q is used in this study, values far less than 10 -20 are acceptable as 
long as an underflow error is not introduced. 

In order to illustrate the effects of the choice of r, one time step of a nonlinear heat conduction problem was 
solved with various choices of r. The choices range from 10 _1 to 10 _lo °, and the results can be seen in Figure 15. 
It is evident that with the larger choices of r, the convergence rate is reduced to first order, but the method still 
reaches the convergence criteria of 10” 12 . The optimal second order nonlinear convergence was reached with values 
of r < 10 5 , and performance remained static with further reductions in r, which is evident from the three co-linear 
results. Additionally, using a factor of 10 _lo ° presented no numerical issues. It is also evident from Figure 15 that 
the iterations required for convergence is reduced from 10 with r = 10 _1 to 6 with r < 10 -5 which is directly 
proportional to the relative computational time for the two cases. From this study, it is evident the the complex-step 
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Number of Elements in One Dimension 


Figure 13. Grid convergence results: MMS gas flow problem. 
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Figure 14. Nonlinear convergence results: MMS gas flow problem. 
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method is a valid way of calculating Jacobians to machine precision which gives the optimal convergence rate of the 
Newton scheme. Unlike calculating Jacobians numerically with other methods such as finite-difference, the choice of 
the perturbation factor for the complex-step method is relatively straightforward. 



Figure 15. Complex-step size variation. 


V. Conclusions and Future Work 

Verification results for the Charring Ablating Thermal Protection Implicit System Solver were presented. Several 
nonlinear terms in the energy and gas continuity equations, as well as many boundary conditions, were verified through 
comparisons with analytical solutions and MMS. Through the verification process, the code exhibited the expected 
second order spatial and second order temporal accuracy. New MMS problems were developed for both the energy 
and gas continuity equations, and their source terms were archived for dissemination to the ablation and heat transfer 
communities. Newton’s method for solving the system of nonlinear equations was verified to exhibit the optimal 
second order convergence for all problems presented, and the complex-step method was introduced to numerically 
calculate exact Jacobian terms. 

Throughout this work, the verification process helped identify several bugs resulting from both coding and deriva- 
tion errors. Consequently, the verification exercises have been invaluable in developing confidence in the new software 
as well as improving robustness, usability, and efficiency. 

Future work will include the development of a MMS problem that exercises the coupling terms in the governing 
equations as well as the decomposition kinetics involved in solving charring ablator problems. No concerted effort 
has been put into improving the efficiency of the software. Consequently, profiling exercises should be performed. 
Validation through comparison with experimental data should be done for the models that have been implemented 
and verified, and the addition of improved physical models, such as ablation models, porous flow laws, and in-depth 
thermochemical non-equilibrium need to be completed. 
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Software Description 



Solves two and three dimensional charring ablation, heat transfer, and porous flow problems with 
moving mesh capability in both serial and parallel 

Uses libMesh finite element library 

•mesh data 

•domain decomposition (METIS) 

•parallel communication (MPI) 

•linear system solver (PETSc) 

Galerkin finite element spatial discretization on unstructured hybrid meshes 
•First order Lagrangian basis functions (second order discretization) 

•Gaussian quadrature 

First and second order implicit time integrators 

Fully implicit linear system 

•Newton’s method 

•Generalized minimal residual method 

Adaptive mesh refinement with all element types 

•2D: quadrilateral and triangle 

•3D: tetrahedron, hexahedron, prism and pyramid 


Uses Cantera for gas thermochemical and transport properties 


Governing Equations 



Implemented in code 


d 


Conduction 


Gas Flow 


Mesh Motion 


dp : 


Solid : —r~ ~ V • {p s v m ) = d> ; 


d 
d{<t>P g ) 


Mesh Motion 


s 

Source 


Gas: — r L + V'(fcv.)-V'(fev,)= co 


d 


Gas Flow 

/C 

Darcy's Law : v ff = 


Mesh Motion 


g 

Source 


<t>P 


VP 



Verified in this study 

tncigy : ■ - V tv 7 + V h \ )- V iph\„ ) = Q 

ut y - v 6 8 , v ' s-< 

Conduction “ M-ct. Source 


Solid :^-V-( A v,„)= d 
a 


Mesh Motion 


Source 


Gas : d ^ Pg ^ + V • ($p g v J- V • {<f>p g v m ) = w. 


dt 


Mesh Motion 


g 

Source 


Darcy's Law : v„ = VP 

<t>ju 
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Linearization 



Linearization in iteration space 



(*>> r), 


+ higher order terms 


Analytically take the derivatives of the residual and implement resulting expressions into 
code 

• Typically approximated for very complicated residual equations 

• Prone to errors in both derivation and implementation 

• Approximation and errors degrade convergence rate of nonlinear system solver 


Finite-difference approximations are also common 

• Eliminates derivation errors 

• Approximations degrade convergence rates of nonlinear system solver 

• Choice of step-size is not straightforward 

• Second order central difference requires 2 perturbed residual evaluations 
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Complex-Step Method 



Taylor series with complex-step 


R 


T. +iAT.. 

j j' 






Real part gives residual 


R 


t p 



= 9i1« 


Tj+ihTj,{fp g )^ 



Imaginary part gives Jacobian 



3- 

[* 

CN 

<1 
* r—i 

+ 

1 1 

fa) 

J} 

AT,. 



Truncation error with extremely small perturbations results in exact residual and Jacobian 
Advantages: Disadvantages: 

• No Jacobian derivation necessary • Extra time spent doing complex arithmetic 

• Easy addition of new models • Incompatible with external libraries 

• Reduces possibility of bugs 

• Reduces development time 

• Optimizes convergence 

• Choice of perturbation size is simple 

• Get residual and Jacobian from one evaluation 
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Verification and Validation 



VERIFICATION : Are we solving the equations correctly? 

• Do the spatial and temporal discretizations exhibit the theoretical orders of accuracy? 

• Does the Newton solver converge at the theoretical rate? 

• Strictly a mathematical exercise 

• Does not need to bear any semblance to the reality of the problems the code intends to 
solve 

• Done through comparison of simulation solutions with known exact solutions 

• Analytic solutions to PDE(s) 

• Method of Manufactured Solutions 

• Should be done for all boundary conditions, governing equation terms, and code options 

• Should precede validation 

• Neither extensive code use nor code-to-code comparisons constitute verification 

VALIDATION: Are we solving the correct equations? 

• Do the model equations represent physical reality? 

• Done through comparison of simulation solution with test data 

• May have to validate and develop best practices for several different regimes/problem 
classes 
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Verification Strategy 



Typically spatial and temporal orders of accuracy are independently verified 

DE = aAt p +bAh q => d[log{DE)] = 

^[log(A/)] 


Simultaneous refinement can be used to verify both temporal and spatial accuracy 


DE = 


1 + 


aAt p ^ 


v bAh q j 

Hold Constant 


bAh q = 



Newton’s method should exhibit quadratic convergence 


E v+l =a{E v } = 


4o g 

u 

^[log(^ v j 





^pf 2-D Energy Equation: Analytic Solution (I) 


Uniform density specimen exposed to 
constant heat flux on one face 

• Adiabatic on three other faces 

• Thermal properties linear in temperature 

• Solution is ID in nature, but suitable for 2D 
and 3D verification with unstructured solver 

Grid refinement study 

• Quadrilateral and triangle meshes 

• 4 grid levels with space/time refinement 

• Coarse grids show here 
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Cg) 2-D Energy Equation: Analytic Solution (II) 




Method of Manufactured Solutions (I) 


Design problem known exact solution 

• exercise multiple terms 

• exercise multiple boundary conditions 

• arbitrary solution domain 

• one solution can be used for multiple verification problems 

• typically non-physical 

Specify PDE(s) to be verified 


P C p — - V-(fcVT) = Q 


Prescribe a solution 



Select property models 
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Method of Manufactured Solutions (II) 

Substitute solution and property models into governing equations and solve for source term 

Q ( x , t) =p [ C po 4 C p j cos(.A x x 4- A t t) cos (B y y 4 B t t) cos (C z z 4 C t t) cos (D t t) 

4 Cpj cos 2 (j4 x x 4 A t t) cos 2 {B y y 4 B t t) cos 2 (C z z 4 C t t) cos 2 (D t t) ] 

{ — 0 [ A t sin(^4 x x 4 A t t) cos (B y y 4 B t t) cos (C z z + C t t) cos (Dtt) 

4 B t cos(a x x 4- A t t) sin {B y y 4 B t t) cos (C z z 4- C t t) cos (D t t) 

4 Ct cos(j4 x x 4 Att ) cos (B y y 4 Btt) sin(C z z 4 Ctt) cos (D t t) 

4 D t cos(j4 x x 4 A t t) cos (B y y 4- B t t) cos (C z z 4 C t t) sin {D t t) ] } 

— { [ k 0 4 fcj cos(A x x 4- A t t) cos (B y y 4 B t t) cos (C z z 4 C t t) cos (D t t) 

4- k 2 cos 2 (j4 x x + A t t) cos 2 (B y y + B t t) cos 2 (C z z + Ctt) cos 2 (Dtt) ] 

[—0 (j4 2 4 By 4 Cl) cos(A x x 4- Att) cos (B y y 4- Btt) cos (C z z 4- Ctt) cos(£>tt)] 

— 9 [ A x sin(.A x x 4 Att) cos (B y y 4- Btt) cos(C z z 4 Ctt) cos (Dtt) 

( — k\A x sin(j4 x x 4- A t t) cos(B y y 4 B t t) cos (C z z + C t t) cos (D t t) 

4 2k? A x sin(j4 x x + Att) cos(A x x 4 A t t) cos 2 {B y y 4 B t t) cos 2 (C z z 4 Ctt) cos 2 (Dtt) ) 

4 B y cos(j4 x x 4 A t t) sin (B y y 4 B t t) cos (C z z 4 Ctt) cos (D t t) 

( — k\B y cos(A x x 4 A t t) sin (B y y 4 B t t) cos (C z z 4 C t t) cos (D t t) 

4 2k 2 B y cos 2 (A x x 4 Att) sin {B y y 4 Btt) cos (B y y 4 Btt) cos 2 (C z z 4 Ctt) cos 2 (Dtt) ) 

4 C z cos(A x x 4 Att) cos {B y y 4 B t t) sin(C 2 z 4 C t t) cos (D t t) 

( — k\C z cos(^4 x x 4 A t t) cos (B y y 4 B t t) sin(C 2 z 4 C t t) cos (D t t) 

4 2k 2 C z cos 2 (v4 x x 4 A t t) cos 2 (B y y 4 B t t) sin(C 2 z 4 C t t) co s(C z z 4 C t t) cos 2 (D t t) ) } 

Implement source term in source code 

Perform verification studies 

•ID, 2D, or 3D 
•Transient or steady state 
•Linear or nonlinear 
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3-D Energy Equation: MMS (I) 




Hexahedro, 


time = 0 seconds 


Tetrahedroi 


time = 20 seconds 


3-D Energy Equation: MMS (II) 




Relative Error at Iteration v 



3-D Gas Flow Equation: MMS (I) 


Gas Flow Equation with constant porosity, permeability, and viscosity 

,dp g kRT 

~ ^ ' (Pg^Pg) ^ g 

Manufactured solution 

p g (x, t) = Ocos(A x x + A t t ) cos (B y y + B t t)cos(C z z + C t t) cos(D t t) + p 9c 


Resulting source term 


uj g ( x , t ) = —4>6 [A t sin (A x x + A t t) cos ( B y y + B t t) cos ( C z z + C t t) cos ( D t t ) 

+B t cos (A x x + A t t) sin ( B y y + B t t) cos (C z z + C t t ) cos (D t t) 

+C t cos ( A x x + A t t) cos (By y + B t t) sin ( C z z + C t t) cos ( D t t ) 

+D t cos ( A x x + A t t) cos (By y + B t t) cos (C z z + C t t) sin (D t t) ] 

t’R-T 

— * — — { — Pg 0 @ ( A x + By + ) cos (A x x + Att) cos (Byy + Btt) cos (C z z + Ctt) cos (Dtt) 

P 

-6 2 ( A 2 x + B 2 + C 2 ) cos 2 (A x x + A t t) cos 2 (B y y + B t t) cos 2 (C z z + C t t.) cos 2 (D t t) 
+6 2 cos 2 (D t t) [A 2 sin 2 (A x x + A t t) cos 2 (B y y + B t t) cos 2 (C z z + C t t) 

+B 2 cos 2 (A x x + A t t) sin 2 (B y y + B t t) cos 2 (C z z + C t t) 





3-D Gas Flow Equation: MMS (II) 



Complex-Step Size Study 



What is the best choice for the step size? AT = rT 

Solved 1 time step of nonlinear conduction problem with 5 choices of perturbation factor 


Large choices still converge 
but at degraded rates 



Very small choices present no 
numerical issues 

Broad range of acceptable 
choices to give optimal 
convergence and reduce 
computational time 

Choice of perturbation factor 
simple compared to finite 
difference 
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Conclusions and Future Work 



Most stationary mesh terms in the software have been verified to exhibit the theoretical 
orders of accuracy 

• Verified both 2D and 3D element types 

• Used analytic solutions and MMS 

Implemented to complex-step method to calculate exact Jacobians 

• Verified convergence rates of nonlinear solver 

• Examined sensitivity to perturbation factor 

Future verification 

• Coupled equation MMS problem 

• Moving mesh terms 

• Adaptive mesh refinement cases 

Modeling additions and software improvements 

• Inverse capability 

• Sensitivity analysis and uncertainty propagation 

• Higher fidelity physical models (thermal and chemical non-equilibrium, etc.) 

• Improved mesh motion algorithms 

• Profiling and optimization 
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